Localised zero-energy modes in the Kitaev model with vacancy-disorder 
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We study the effects of vacancy disorder on the Kitaev model defined on a hexagonal lattice. We 
show that the vacancy disorder induces a zero-mode that is localized at the defect site. We derive 
analytical forms for these localized wave functions in both the gapped and gapless phases of the 
Kitaev model. We conjecture that the vacancy disorder can be utilized as a probe of the quantum 
phase transition (from the gapped to gapless phases) in this model. The behavior of the Inverse 
Participation Ratio (IRR) in the gapless phase and across the transition is also studied numerically. 
Comments are made about the behavior of site-site entanglement in the single particle states for 
the case of a single vacancy. 



I. INTRODUCTION 

The effect of quenched disorder typified by impuri- 
ties, lattice imperfections, and vacancies on condensed 
matter systems has been a source of intense scientific 
investigation in the recent past. In fact, the study of 
such "frozen-in" disorder has led to the unravelling of 
a host of very interesting phenomena like infinite ran- 
domness fixed points [H , quantum Griffiths effects [l|, [|J , 
and maybe even smearing of phase transition Even 
though as highlighted above multi-impurity effects can be 
extremely interesting, the study of a single impurity em- 
bedded in a host can also act as an efficient probe of the 
physical characteristics of the underlying bulk material, a 
situation exemplified by the case of local impurity acting 
as a probe of the order-paramatcr symmetry in uncon- 
ventional superconductors This manuscript for the 
most part belongs to the latter genre wherein we study 
the role played by a single impurity in identifying the 
quantum phase transition inherent in the Kitaev model. 

The Kitaev model has become one of the paradigmatic 
models that has been studied in various contexts ranging 
from strong correlation physics to topological quantum 
computation. Its theoretical appeal lies in the fact that 
it represents one of the few spin systems that can be 
solved exactly. The solution to the clean Kitaev model is 
effected by recasting the spin Hamiltonian into that of an 
equivalent Majorana hopping problem in the background 
of static Zi gauge configurations. The exact solution of 
the model reveals both a gapless and gapped spin-liquid 
phases with a zero-temperature quantum transition in- 
terpolating between these two phases. The gapless spin 
liquid phase is quite unique as it supports a spin-spin 
correlation function that is short ranged [6j , thus setting 
it apart from other spin-liquid phases studied so far. It 
also supports fractionally charged topological excitations 
both Abelian and non-Abelian that can be plausibly used 
to perform quantum computations. Apart from its util- 
ity for topological quantum computation or for its useful- 
ness in studying spin-liquid ground states, the spin-1/2 
Kitaev model defined on a two dimensional hexagonal 
lattice has become a powerful test-bed example to study 
various fundamental concepts in the field of strong cor- 
relation physics. For instance, it has been used to study 



fractionally charged excitations that occur in topological 
insulators thereby providing a beautiful higher dimen- 
sional extension Q of the Jackiw-Rebbi theory [l2| , [H[ > 
that describes charge fractionalization in one dimensional 
system. It has also been utilized to study dynamics of 
quantum quenches across the critical region Q- More- 
over, there now exists higher dimensional realizations of 
the Kitaev model [l6[ , and also extension to higher values 
of spin [l7[ 

As is clear from the preceding paragraph the clean Ki- 
taev model has been the subject of intense scientific in- 
vestigations in the last few years. However, apart from a 
few notable exceptions, (detailed in the next paragraph), 
one area that has been rather neglected is the study of 
the effect of impurities on the Kitaev model. This is a 
particularly glaring deficiency as now there exists pro- 
posals for experimental realizations of the Kitaev model 
[Tol | . [Hj], @- Thus, the study of impurity effects gain an 
added significance as realistic systems are seldom clean. 

Now, the simplest form of quenched disorder involves 
studying the effect of a single impurity on the bulk sys- 
tem. Such a study was undertaken by [i"oT | for the case 
of a single magnetic impurity that was embedded in the 
host Kitaev model. It was shown that coupling of an 
impurity to the host Kitaev system leads to an unusual 
Kondo effect that is sensitive to the topological transition 
in the Kitaev model. In a related work, Willans et. al 
[l4| showed that disorder in the form of a single vacancy 
binds a flux which in turn gives rise to a local moment. 
Furthermore they showed that this moment leads to a 
vacancy susceptibility that diverges logarithmically as a 
function of the applied field, (for weak applied fields). 

While it is true that impact of quenched disorder has 
received very scant attention, the effect of impurities on 
allied models have been rather well studied. More specif- 
ically: The Kitaev model maps onto a fermionic model 
displaying bi-partite hopping on a hexagonal lattice in 
the background of Z2 gauge fields. Now, the impact of 
quenched disorder on similar bi-partite hopping models 
have been studied in the context of Anderson localiza- 
tion. In a seminal work, (see 18 1 [Hj]), it was shown that 



the quenched impurities lead to a divergent Density of 
States, (DoS) in models that display bi-partite hopping. 
More specifically, by using a field theoretical formulation, 
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these authors showed that a random-mass form of dis- 
order, (in addition to a random vector potential) would 
lead to a highly divergent DoS that conforms to the func- 
tional form, p(E) ~ -ie-l 1 "^ 17 *, with x = 2. This faster 
than power-law divergence of the DoS should be con- 
trasted with the results of Ludwig et. al [13] wherein they 
studied a random vector potential model with bi-partite 
hopping. These authors showed that the DoS in these 
models diverge as a power- law p(E) ~ E~ 1+2 / z , where z 
is a continuously varying dynamical exponent. Now, the 
field theoretical treatment of Gade and Wegner suffered 
from a slight draw-back: It did not provide a physical 
frame-work wherein the diverging DoS at the band cen- 
ter could be understood. This situation was remedied in 
the work of Motrunich et. al (2lj j wherein an intuitive 
and physically appealing argument was provided for the 
origin of the diverging DoS and the low-lying states that 
was causing it. They further argued that the DoS in- 
deed does diverge with a functional form analogous to 
that derived by Gade and Wegner, however, with the ex- 
ponent x — 3/2 instead of 2 obtained in Another 
extremely relevant physical context wherein such disor- 
dered bi-partite hopping models holds relevance is pro- 
vided by the case of graphene. In graphene, due to open 
surfaces and substrates, disorder is an ever present bug- 
bear. One form of disorder that has been relatively well 
studied in the context of graphene is vacancy disorder. 
This type of disorder arise naturally in the case of irra- 
diated graphene wherein the carbon atoms are knocked 
out of graphene planes. The impact of such vacancies on 
the electronic properties of g raphene were investigated 
in a series of papers H|, Haj. The influence of various 
other forms of disorder, (also inclusive of the case of va- 
cancy disorder) on the electronic properties of graphene 
was studied in [22j], [23j |. 

In this paper we will study the influence of vacancy 
disorder in the Kitaev model. More specifically, we focus 
our attention mainly on the structure of the wavefunc- 
tion at the site of a vacancy disorder. We show that 
vacancy disorder gives rise to a "zero-mode" that is lo- 
calized at the impurity site. As we shall see later on in 
this paper, this zero-mode exists in both the gapped and 
gapless phases of this model and is a consequence of the 
particle-hole symmetry of the bi-partite hopping prob- 
lem. However, as will be shown in the bulk of this man- 
suscript, the functional form of the zero- mode is quite 
different in two phases thereby providing an invaluable 
tool for distinguishing different phases of this model. The 
interesting question of impact of many vacancies on the 
Kitaev model will be addressed in a future publication 
24|. 

The paper is organized in the following manner: In 
Sec. |TT] we recapitulate the mapping of the Kitaev 
model into a non-interacting Majorana fermion problem 
by following the Jordan- Wigner fermionization scheme. 
Sec. IIIII is devoted to the analytic derivation of the zero- 
mode wave function that is localized at the vacancy site. 
The functional form for the wave function is established 




FIG. 1. Hexagonal lattice with a vacancy. A unit cell contains 
two points, one point marked red and another unmarked point 
forming a bond. The unit cells are labelled using pair of 
integers (j, I) and the convention used is made clear through 
explicit labeling of some unit cells. The vacancy is denoted 
by a yellow circle. The dotted line separate the part of the 
lattice with j > from the part with j < 0. 



for both the gapless phase. UlI Al and for the gapped spin- 
liquid phases ITlI Bl In MI CI we will delve into the issue 
of the "flow" of the Inverse Participation Ratio (IPR) of 
the localized wavefunction as a function of the coupling 
parameters in the Kitaev model. In the same subsection 
we will also touch upon related entanglement measures. 
Finally, we will end with a concluding section, Sec. [V] 
wherein our results will be briefly reviewed and placed in 
context of the existing literature in this field. We will also 
briefly mention some of the open problems that remain 
to be tackled in this subsection. 



II. KITAEV MODEL 

In this section we give a short introduction to the 
model and its properties for later reference. It also serves 
to set the notations for the rest of this paper. The model 
comprises of spins residing on the sites of a honeycomb 
lattice as shown in Fig. [T] The spins interact with each 
other via nearest neighbour coupling which is dependent 
on the bond orientation. These orientations are labeled 
as as x, y and z in the figure, Fig.[TJ Also, as represented 
in Fig. [1] the y-link is taken to be the basis, with a two 
"atom" unit cell: The red colored lattice point denot- 
ing the A-sublattice and the uncolored point indicative 
of the B-sublattice respectively. They are connected to 
each other via the y-link. A point in the lattice is thus 
labeled by a triplet of numbers (j, I, fx) where j, / denote 
the unit cell and p = 1(2) correspond to the A(B) sublat- 
tice. Thus, under this labeling scheme the Hamiltonian 
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is expressed as: 

+ Jz<Tj,l,l<Tj-l,l-l,2\- (!) 

Here, as usual the as are the usual Pauli matrices that 
represent the spin variables. 

As briefly discussed in the introduction, Sec. U this 
spin-model can be mapped onto a Majorana Fermion 
hopping problem. Different methods can be adopted to 
effect this transformation. In this manuscript we shall 
employ the Jordan- Wigner fermionization scheme as em- 
ployed by Feng et. al., 26] in this context. Define the 
Jordan- Wigner tail operator as 

K(j,l,fl)= [] <n,», (2) 

( t 7,/,/i)>(m,n,^) 

where (j,l,fj,) > (m, n, v) if j > m or j = m,l > n 
or j — m,l — n, /i > v. Now the Majorana fermion 
operators can be defined as: 

1>h, M = K(j, I, li)of tliM , = K(j, I, »)al^. (3) 

In terms of these operators the Hamiltonian takes the 
form 

3,1 

+J z D^l L1 ^_ hl _ h2 \. (4) 

Here, the operators Dji = itpj ; j^-i i_i 2 , defined on the 
•z-links, is Hermitian and commutes among themselves 
and with the Hamiltonian reflecting the local symmetry 
of the Kitaev model. It can be shown that the operators 
Dji have eigenvalues ±1. The Hamiltonian gets block 
diagonalised into different sectors corresponding to dif- 
ferent sets of eigenvalues of Dji. In each of these sectors, 
the Hamiltonian becomes a quadratic fermionic system 
obtained by replacing each Dji by its eigenvalue and can 
be re-cast into the form: 

H = ^ T iA^, (5) 

where if} = (■•' V'j i iVs i 2' •■) T an d A is an antisymmetric 
matrix. 

Thus, as alluded to in the introduction the Kitaev 
model has been mapped onto a non-interacting Majorana 
fermion problem in the background of static Z 2 gauge 
field. 

If there are N number of unit cells, we have 2N spins 
and the Hilbert space is 2 2N dimensional. Since there 
are N z-links, there are 2 N sectors each of which has 
dimension 2 N corresponding to 2N Majorana fermions. 

For further calculations, let us first see how the eigen- 
vectors/values of the coefficient matrix iA are related to 
the fermionic excitation modes of the system. Since A 
is antisymmetric, the eigenvalues of iA comes in pairs 



— ej, €i with eigenvectors Vi, v* respectively, where ti > 0. 
We can choose the eigenvectors to be orthonormal since 
iA is Hermitian. Define fermion operators di — ^^Vi- 

It is easily seen that these operators obey {d\,dj} — Sij. 
We get, 

N f 1\ 
# = \Adi-~ . (6) 

It is known that the ground state of the Hamiltonian 
lies in a sector wherein all the Dji operators take the 
eigenvalue +1 [2|| |27|- By making use of the transla- 
tional symmetry in the problem, a solution of the model 
can be effected by going into Fourier transformed rep- 
resentation . Thus, the Hamiltonian re-expressed in 
terms of the Fourier transformed variables, (V'k V'k) = 
E A1 e- lk ' rj ' ; 1#,i, 2 )/V2JV, (where r jJL = jm + ln 2 ) 

reads, 

» -\ Efc A) (^ViT )(!)■<') 

Here, </>(k) = 2(J x e- lk2 + J y + J z e'^ kl+k ^) and h = 
k • n^. The eigenvalues are ±|0(k)| and the fermionic 
excitations are given by 

# =5>(k)i(dt(kMk) iy (8) 

where rf(k) = ^[V'k + i |0(k)| V'k]- ^he excitation spec- 
trum is gapless if there exist points where |<^(k)| = 
which is possible only if following condition is satisfied: 

{J X ~ Jyf < Jl < (J X + Jy) 2 . (9) 

The gapless phase is characterised by Fermi points where 
the 4>(k) vanishes. There is a quantum phase transition 
from the gapless- to the gapped-phase as the parameters 
cross the conditions in Eq.[9l This quantum phase transi- 
tion is the one that we wish to probe via a single vacancy 
disorder. 



III. SINGLE VACANCY 

In this section, we study the nature of the wave- 
function at the site of a single vacancy. The analytic 
functional form is derived for both the gapped and gap- 
less phases. 

For the sake of concreteness consider the Kitaev model 
with a vacancy at the B-site in the unit-cell (—1,-1), (see 
Fig. [T]) . Note that the Jordan- Wigner construction goes 
through with the tail operator K(j, /, /z) missing crf.j _ 1 2 
for all (j, I, /i) > (—1, —1, 2). As the system is no longer 
translationally invariant one cannot use Fourier trans- 
form to solve the problem. However, the general struc- 
ture of the Hamiltonian, Eq.[U remains with A now being 
a (27V— 1) x (2N— 1) matrix obtained by removing the row 
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and column corresponding to the site (—1,-1,2) from the 
matrix A in Eq.[5] Thus, we have TV— 1 eigenvectors form- 
ing pairs as described in Sec. |TT] and one unpaired eigen- 
vector denoted by v. This eigenvector should be real with 
zero eigenvalue because of the e •<-> — e symmetry briefly 
alluded to in SecHH The N — 1 pairs can be combined to 
form N — 1 complex fermion operators di leaving behind 
a single unpaired mode. This unpaired eigenvector forms 
a Majorana mode, d = i\) T v as <v = d and d 2 = 1. Note 
that by removing a spin at (— f , — f , 2), we have left out 
one other Majorana fermion operator from the Hamilto- 
nian; V'o o l wmcn would have formed the operator Dq^ 
with ip°Li _i 2 had the spin been present. Thus a com- 
plex fermion mode can be constructed from these two 
Majoran modes as dN = ^(d + iipQ i) which will be a 
zero-energy excitation of the Hamiltonian which, again, 
has the same form as Eq. |6l However, since Dqq is not 
present in the Hamiltonian, the number of Dji operators 
is now N — 1. Thus, now there are 2 JV_1 sectors each 
with N complex fermions. Hence as expected the total 
dimension of the system is 2 2N ~ 1 . 

Now that the above discussion has clearly established 
that a single vacancy induces a zero-energy fermionic ex- 
citation mode in the Kitaev model, let us turn our at- 
tention to the analytic structure of these modes. To do 
so we employ a method developed by Pereira et. al [28j . 
[29j in the context of zero modes arising out of a vacancy 
defect in graphene. This adaptability of the technique de- 
veloped for the case of graphene to the Kitaev model is 
not so surprising as they both give rise to similar fermion 
hopping problems. Unlike the case of graphene studies 
by Pereira et. al. [HI], j^, where one is restricted to the 
isotropic case J x — J y — J z , here we consider the general 
anisotropic hopping problem and obtain expressions for 
the zero mode for the parameter regimes corresponding 
to both gapped- and gapless-phases of the clean model. 
More specifically, we obtain an asymptotic form for the 
defect wave-function in the gapless phase, whereas one 
can evaluate an exact form of the wave function in the 
gapped phase. 

Before obtaining explicit expression for the zero modes 
of iA with B-site vacancy, let us first see how they are re- 
lated to the corresponding zero mode when the vacancy 
site is in the A-sublattice. Let us introduce the nota- 
tion r ee (j r ,l r ) and A Atj „(r,r / ) ee A^^^^,^,^. The 
clean model has the symmetry given by t x A(y + p, r + 
p')r x — —A{v — p, r — p') for any r, where t x is the 
Pauli matrix. But the vacancy breaks the translational 
invariance. Let V^ ro be the matrix to be added to A 
to create the vacancy, by removing corresponding matrix 
elements from A, at position ro in the sublattice /i. Now, 
'r'Vi.rofro + p,r + p ')t x = -V 2 , ro (r - p,r - p'). If 
^2, r ( r ) = </ ) ( r — r o) is an eigenvector of i(A + V2, ro ) with 
eigenvalue A, then it follows that 0i, ro (r) ee r x (f)(r — r) is 
an eigenvector of i(A + VI, ro ) with eigenvalue —A. Thus 
we need to find the localised zero mode with one type of 
vacancy only, the other obtained from the relation given 



the zero mode created by /i-sublattice vacancy at site r. 

Consider a B-site vacancy in the unit cell (—1, —1) as 
shown in the Fig. Q] The eigenvalue equation of iA for 
the zero eigenvalue decouples the A- and B-sublattice 
amplitudes. Denoting the A-sublattice amplitude by a,ji, 
we get the corresponding eigenvalue equation as 



Jydjl + J x ajJ + l + J z Qj + lA + l — 0. 



(10) 



This equation hold true everywhere except for j = — 1 = 
I , where it no longer applies due to the vacancy. A similar 
equation can also be written down for the B-sublattice 
amplitudes. In that case, the corresponding equation is 
satisfied by choice of them being equal to zero identically. 
Thus, we have ^2,(-i.-i)( r ) = {aj r i r 0) T . 

To solve Eq. [TUl following the procedure of Pereira et. 
al. [29j], the lattice is divided into two parts: j > 0, and 
j < 0, the parts that lies below and above the dotted 
line respectively in Fig. [T] Eq. [TU] is solved separately 
in these two regions and a boundary matching condition 
is imposed at the dotted line in Fig. [TJ Applying peri- 
odic boundary condition along the horizontal direction, 
a Fourier transformation , aj(q) — ^2ie~ zql o,ji, reduces 
Eq.Hnito 



The solutions are given by 



dj(q) 



{-f(q)] j a (q) V j > 0, 
-f(q)Y + 1 a^(q) Vj<-1. 



(11) 



(12) 



Seeking solutions that decay as a function of the distance 
from the vacancy site, we get the following conditions: 
ao(q) is non-zero only if \f(q)\ < 1, and a_i(q) is nonzero 
only if | f(q) \ > 1. Note that these conditions require that 
the j < and j > regions have complementary sets 
of wave-vectors contributing to the eigenvectors. The 
boundary condition at the interface is now implemented 
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ao(q) 



7. 



a-i(q) 



0, (13) 



above. In the discussion hereafter, 



always denote 



except for I = —1. This set of equations is satisfied by 
the choice a (q) = 6(1 - |/(g)|) and (J - + y~ W) 0-i(g) = 
O (| f(q) | — 1). It is easily checked that the condition 
1/(9)1 < 1 can be satisfied for parameter values that obey 
Eq. [9] For other parameter values, corresponding to the 
gapped-phase of the clean model, |/(<7)| will either be less 
than 1 or will be greater than 1 for all values of q, thereby 
giving us trivial solution in one of the two regions. We 
now consider these two parameter regions separately. 



A. Gapless phase 

First, consider the gapless phase, where we have a 
set of q values in the range (q* ,2ir — q*) that satisfy 
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1/(9)1 < 1) where cos(<f ) 



7 — 7 — 7 

2 J x 3 it 



with q* £ (0,7r). 



Its complement in [0, 2ir) gives the set of q values con- 
tributing to eigenvector in the j < region. The eigen- 
vector for the j > region is now constructed by taking 
the inverse Fourier transform of aj(q). Thus, we have 

a 3l ~ K {_£ dqeW-JM [-e(q)e ie ^] j } . (14) 



Here, e(q) = [J x 
tan(%)) = (J x 
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2J x J y cos(q)} 1 / 2 /J z and 
J y )t&n(q/2)/(J x + J y ). Notice that 
e(g) decreases monotonically from its maximum value 1 
at q* to \J X — J y \/J z at q — n. For asymptotically large 
values of j, the dominant contribution comes from the 
region around q* . Therefore, expanding around q*, the 
above equation, Eq. Q3] can be written in terms of its 
asymptotic form as 



a(x, y) ~ di ■ 



e i 9 *3:/v / 3+j2(7r+e*)i;/3 

ay j V3 — i2/3j//v3 — ia; 



(15) 



The parameters a, 0*, and /3 are given by a — 
2J x J y sm(q*)/Jl 9* = %*),/? = (J 2 - J 2 )/J 2 . Also, in 
the above equation x(j, I) — y/3(l — j/2) and y(j, I) = |j 
are the re-defined lattice indices. The integral in Eq. [T4l 
vanishes as we approach the boundary, J 2 — > {J x — J y ) 2 , 
of the parameter regime defining the gapless phase, since 
q* — » 7r here. In the opposite limit of J 2 — > (J x + J y ) 2 , 
q* — > 0, and as J z crosses this condition we move into 
the gapped phase solution which will be discussed in the 
next section, Sec. IIIIBI 

For J x = J y = J z , we have q* = 2ir/3, a = y/3, 9{q) = 
0, and thus the result of Pereira et.al., 29], a(x,y) ~ 
^^2^x/3V3+i2wy/3y^ y _ ix ^ is recoverec i. A numer- 
ically exact zero mode in the gapless phase is shown in 
Fig. @ for a finite system with periodic boundary con- 
ditions. 



B. Gapped phase 

Now consider the gapped phase. For the sake of 
concreteness, consider the situation wherein all Js are 
taken to be positive and furthermore satisfy the condi- 
tion J x + J y < J z . As we have already seen, we have 
| /(g) | < 1 for all values of q and hence the solution for 
j < is trivially zero. Then the boundary condition at 
j = implies (see Eq. [TUJ) that 

aai = V I ^ 0. (16) 

The solution of Eq.[TU] satisfying this boundary condition 



a 0,0i 



(17) 



for all j > 0,1 £ {0, .., j} and zero everywhere else. Here 
iCi is the Binomial coefficient. We note that the solution 




FIG. 2. The zero mode intensity |a.,z| 2 for the gapless case 
(Jx = Jy = Jz = 1) (top) and for the gapped case (J x — J y = 
1, J z = 2.05) (bottom) for a system with iV = 2500 unit cells 
and with periodic boundary conditions. 



is non-zero only in a cone-shaped region extending in the 
j > direction. For any j > 0, 



Oo,0- 



(18) 



Therefore \dji\ decay exponentially since (J x + J y )/ J z < 
1. Note that we have implicitly assumed that the lattice 
extends infinitely in the j direction. If we have periodic 
boundary condition in the j direction as well, the tail of 
this solution can wrap around to the j < region shown 
in Fig. [TJ The zero mode corresponding to the gapped 
phase, but still close to the transition (at J z = J x + J y ), 
is shown in Fig. (|2|). 

We have two other possibilites, namely, J y + J z < J x 
and J z + J x < J y , for the gapped phase. They also 
give similar results, and are related to the current result 
by rotation of the lattice by 2ir/3 and 4ir/3 and cyclic 
permutation of J x ,J y ,J z - 



C. Participation ratio and site-entanglement 

The contrasting nature of the zero modes in the gapped 
and gapless phases provides motivation for a closer study. 
In the gapless phase there is a "quasilocalized" zero mode 
in the terminology of [28| . as the amplitude decreases as 
1 jr from the vacancy. This leads to an anomalous scaling 
of the inverse participation ratio (IPR) defined as 



P 



E 



(19) 
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The IPR in the gapless phase would then depend on the 
size of the system N as l/ln(iV) 2 [28|, whereas in the 
gapped phase the IPR would be independent of the sys- 
tem size reflecting the localized nature of the zero mode. 
In Fig. ([3]) is shown the IPR across a transition to the 
gapped phase where we can see an increased localization, 
as indicated by the rapid increase in the IPR beyond 
J z = 2. 

Quite apart from this dependence, it is interesting to 
see strong variations of the IPR within the gapless phase 
as a function of the parameters (J x , J y , J z ). In Fig. ([3]), 
top panel, this is seen in the region J z < 2. Also note 
the strong dependence of these oscillations on the system 
size N in this case. The variation of the IPR in the en- 
tire gapless phase is most neatly captured in the triangle 
with J x + J y + J z = 1, with all J x , J y and J z being < 
1/2 [25j]. The IPR of the zero mode for parameter val- 
ues in this triangle corresponding to the gapless phase 
is shown in Fig. ([3]). The borders of the triangle corre- 
sponding to an imminent transition to the gapped phase 
shows a minimum of the IPR, indicating the existence of 
extended zero modes. The dark regions corresponding to 
very small IPR and large derealization are arranged in 
an intriguing manner and require further work for elu- 
cidation. The complexity of the figure in terms of the 
number of such regions with large derealization increases 
with the system size N. 

A trivial calculation on the gapped side shows that the 
wavefunction is not only square summable but actually 



summable: 



< 00. As the gapless phase bound- 



ary is reached this summability is lost. The total site- 
entanglement present in the one-particle modes is closely 
related to this sum. Entanglement in the Kitaev model 
has been recently studied, and refers to entanglement in 
the spins [30] . However if we are to look at single parti- 
cle states, we can study entanglement between the sites 
themselves, sites that maybe empty or singly occupied, 
and the mode is considered to be a superposition of such 
singly occupied states. 

In the context of the Kitaev model the onsite fermions 
are of Majorana type as opposed to complex ones. Al- 
though the site-entanglement measure in the context of 
Majorana fermions needs further studies, on interpreting 
the modes as that of a complex fermion hopping prob- 
lem, site entanglement becomes an especially standard 
and well studied tool. Such entanglement measures have 
been used previously in many contexts including that 
of Anderson localization [3lj |. wherein a single site von 
Neumann entropy has been studied. If however we study 
the entanglement between a pair of sites, say labelled 
by (j, I) and (f , I'), the concurrence [32[ measure can be 
used. The concurrence measures entanglement between 
any two two-state (qubits) systems, and sites with occu- 
pancy or 1 are precisely isomorphic to qubits. If the 
concurrence is 0, there is no entanglement between the 
two sites and if it is 1, they are maximally entangled. For 
one-particle states the concurrence is simply 33] 




FIG. 3. 



The inverse participation ratio of the zero modes 
when J x = J y = 1 as a function of J z (top), the gapless to 
gapped transition being at J z = 2. The IPR as a function 
of J x , J y and J z in the entire gapless phase (bottom) for 
N = 900. The darker (blue color) regions have a low IPR or 
large participation ratio. 



and the total concurrence, summed over all pairs of sites 



is Ct where 



Ct — 




■1 = a 



J : 



0.0 



{J z Jy Jx) 



-1. (21) 



Cji.j'v = 2|aj ) ja J v ji /| 



(20) 



Thus, at the transition when J z = J x + J y we see a di- 
vergence of the total site-site entanglement. The inverse 
participation ratio P is also simply related to site-site en- 
tanglement. The sum of the squares of the concurrence 
(also called the tangle) across all pairs of sites is related 
to the IPR. While a closed form analytical expression 
for the IPR seems difficult, as noted above when dis- 
cussing Fig. ([3]) the IPR is a minimum across the gapless- 
gapped transition, indicating increased derealization of 
the states and large site-site entanglement. 



IV. TWO VACANCIES 

Next we briefly discuss the effects of having a vacancy 
pair. To do so, let us first consider a sort of index theo- 
rem given in Pereira et. al., in the context of graphene 
p8j . Generally, this "index"- theorem counts the num- 
ber of zero modes that arise due to presence of vacancy 
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defects in a fermion hopping problem on a bi-partite 
lattice. More specifically, it has been shown that the 
number of vacancy induced zero modes in such tight- 
binding type models is equal to the difference \ub — tia\, 
wherein 713(71,4) is the number of vacancies on the A(B) 
sub-lattice. Now, it is also known that for instance if 
riB > 77^, then the zero modes have non-zero support 
only on the A-sublattice. The situation is reversed if 
riB < tia- Thus, according to the above discussion if the 
vacancy pair is introduced on different sublattices then 
one would assume that the zero-modes interact with each 
other lifting away from zero. While this is indeed true 
in the gapless case the gapped case comes with an addi- 
tional wrinkle. In other words, in an infinite lattice with 
open boundary condition, depending on the position of 
vacancies, there may still be intact zero modes even when 
the two impurities are placed on different sub-lattices. To 
see this let us define A to t — A + V\, ri + V2, r2 . Then 

{A tot 4> 2 ,T 2 ) (r) = 2</> a (ri 2 )[J !; <5 r ,r 1 + ^x^r,n-(o,i) 



-«7*<Sr,ri-(l,l)] I l 







(22) 



where we have written <\> = (<fr a 0) T and 1*12 = ri — r 2 . 
Comparison with the previous section gives <p a ( r ) — 
a j r -i.i r -i- For the gapped phase, 4> a (^i2) is zero unless 
ri is within the cone where the zero mode is nonzero. 
Thus, 02, r 2 is also a zero mode when ri is outside this 
cone. Note that the fact that the lattice is infinite in 
the j direction is crucial for this argument. For peri- 
odic boundary condition, the cone could wrap around 
the torus and the position ri will be within the cone. 
For other cases including the gapless phase, a (i"i2) is 
nonzero in general and we could represent the effective 
coefficient matrix in the space spanned by the two zero 
modes 



iA 



iS 
-iS 



(23) 



where S = 0i )ri Vi, ri 02,r 2 = 2^ o (r 12 )[J J/ ^ o (0) + 
Jx<l>a((0, 1)) + J z 4> a {{l, 1))]. Here we have also used the 
relation <\>\ ' Vi, ri 02,r 2 — — ^I^^.^^i.ri which follows 
from </>i. ri (r) = t x 4>(ti — r) and the relation between V\ 
and V2. The quantity within the square brackets is ex- 
actly what is excluded from being zero in Eq. [TO] S is real 
since the zero modes are real. The eigenvalues of iA are 
±S and the two zero modes lift off from zero eigenvalue 
and a crude estimate of the new eigenvalues of iA to t is 
given by ±5(ri 2 ) that decays with the distance between 
the two impurities: S decays as powerlaw, asymptoti- 
cally, with ria in the gapless phase and exponentially 
when the A-site impurity is inside the cone defining zero- 
mode eigenvector for B-site impurity. 



V. CONCLUSION AND OPEN PROBLEMS 

The role played by vacancies in identifying the gapped 
and gapless phases has been discussed. In particular it 



has been shown that a single vacancy in the gapless phase 
leads to a "quasi-localized" zero-mode that asymptoti- 
cally decays as a power-law. In the gapped phase, the 
zero-mode due to the vacancy defect is exponentially lo- 
calized with-in a cone that emanates from the vacancy. 
These results were obtained analytically by laying re- 
course to a technique developed in the context of [28| . 
[29l | . This leads us to conjecture that a single vacancy im- 
purity can act as a probe in distinguishing the two-phases 
of the Kitaev model. These two phases are characterized 
by very different behaviors of the IPR as well, and while 
the transition is characterized by a local minimum of this 
quantity, it shows for finite lattices, intriguing patterns 
as a function of the parameters in the gapless phase. The 
localization in the gapped phase leads to summable wave- 
functions and to a finite total site-site entanglement as 
measured by the concurrence. This diverges as the gap- 
less phase is approached in a manner that is very easy to 
calculate. 

We have also briefly discussed the effect of interacting 
zero modes. More specifically, specializing to the case 
of two vacancy defects, we have seen that the number 
of zero-modes in the gapless phase is equal to the dif- 
ference of vacancies in the A and B sub-lattice, in con- 
formity with the "index" -theorem in [28], [29]. We have 
also argued that there are situations in the gapped phase 
of the Kitaev model, (in the infinite lattice limit with 
open boundary conditions), wherein the above mentioned 
"index" -theorem does not hold. 

Now, we turn our attention to some open problems 
that still remain to be addressed with regards to the 
effect of vacancy disorder on the Kitaev model. As is 
obvious from this paper, a proper discussion of multi- 
vacancy effect in the Kitaev model is sorely lacking. As 
a prelude to any such effort one needs to generalize the 
zero-mode counting argument of Pereira et. al. so as 
to account for the extreme directionality dependence of 
the zero-mode wave function in the gapped phase. In 
the limit of multiple impurities it is plausible that one 
can reduce the problem to a system that is governed by 
an effective frcc-fcrmion action with both random mass- 
term and random Z2 gauge fields. It is an open question 
whether this is indeed the case. If one could write down 
such an effective Hamiltonian, in the spirit of [2(3], fl9j . 
Jl8], one could analytically investigate the effect of impu- 
rities in determining thermodynamic properties like the 
DoS. It would be interesting to see whether these models 
show Griffiths type behavior exemplified by a divergent 
DoS wherein the divergence is controlled by continuously 
varying exponents that are a function of the disorder con- 
centration, (see Ref. 34{ where a similar effect was shown 
to exist in the ±J random bond Ising model). Results 
that come from such effective action description of dis- 
order effects can also serve to shed light on the effect of 
disorder on spin-liquids in general 

Some of these issues addressed above could be also 
studied numerically. More specifically, the functional 
form of the DoS as a function of the disorder concen- 
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tration and other system parameters are being studied 
by numerical investigations 24] . 

As this paper was being written up, we were made 
aware of a pre-print [35j wherein results similar to ours 
in the context of Kitaev model were obtained. 
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